La base de données utilisée dans ce TP concerne la géolocalisation des accidents de la route à Paris. Il s’agit plus précisément des bases de données annuelles des accidents corporels de la circulation routière, en particulier le millésime 2019.
« Pour chaque accident corporel (soit un accident survenu sur une voie ouverte à la circulation publique, impliquant au moins un véhicule et ayant fait au moins une victime ayant nécessité des soins), des saisies d’information décrivant l’accident sont effectuées par l’unité des forces de l’ordre (police, gendarmerie, etc.) qui est intervenue sur le lieu de l’accident. Ces saisies sont rassemblées dans une fiche intitulée bulletin d’analyse des accidents corporels. L’ensemble de ces fiches constitue le fichier national des accidents corporels de la circulation dit « Fichier BAAC » administré par l’Observatoire national interministériel de la sécurité routière “ONISR”.
Un certain nombre d’indicateurs issus de cette base font l’objet d’une labellisation par l’autorité de la statistique publique (arrêté du 27 novembre 2019). »
Vous pouvez télécharger les bases de données brutes ici ou utiliser R pour les télécharger dans votre dossier actuel :
# download the dataset
download.file("https://github.com/comeetie/quantilille/blob/master/exercises/data.zip?raw=true",
destfile = "data.zip")
# unzip
unzip("data.zip",exdir=".") sf et les data.frame associéssf::st_read().
library(sf)
iris.75 <- st_read("data/iris_75.shp")Reading layer `iris_75' from data source `C:\Users\Kim Antunez\Documents\Projets_R\quantilille\exercises\data\iris_75.shp' using driver `ESRI Shapefile'
Simple feature collection with 992 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 643075.6 ymin: 6857477 xmax: 661086.2 ymax: 6867081
Projected CRS: RGF93 / Lambert-93
plot(iris.75). Que remarquez-vous ?
plot(iris.75)sf.
sf::st_geometry() ? Quelle solution au problème précédent proposez-vous ?
sf::st_geometry() permet d’isoler l’information contenue dans la colonne geometry de l’objet sf. Cela permet de mettre de côté les autres variables et de n’en afficher qu’une.
plot(st_geometry(iris.75))plot.
sf::st_read() et sf::st_geometry(). Vous pouvez aussi customiser la carte en utilisant différents paramètres de la fonction plot : bg, col, lwd, border, pch, cex…
accidents.2019.paris <- st_read("data/accidents2019_paris.shp")Reading layer `accidents2019_paris' from data source `C:\Users\Kim Antunez\Documents\Projets_R\quantilille\exercises\data\accidents2019_paris.shp' using driver `ESRI Shapefile'
Simple feature collection with 11897 features and 7 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 643222.2 ymin: 6857482 xmax: 661451.8 ymax: 6867053
Projected CRS: RGF93 / Lambert-93
plot(st_geometry(iris.75), bg = "cornsilk", col = "lightblue",
border = "white", lwd = .5)
plot(st_geometry(accidents.2019.paris), col = "red", pch = 20, cex = .2, add=TRUE)
title("Accidents à Paris")Utilisez sf::st_intersects() et sapply().
Documentation de la variable grav : Gravité de blessure de l’usager, les usagers accidentés sont classés en trois catégories de victimes plus les indemnes :
1 : Indemne 2 : Tué 3 : Blessé hospitalisé 4 : Blessé léger
library(dplyr)
inter <- st_intersects(x = iris.75, y = accidents.2019.paris)
inter_nonbless <- st_intersects(x = iris.75, y = accidents.2019.paris %>% filter(grav==1))
iris.75$nbacc <- sapply(inter, length)
iris.75$nbaccnb <- sapply(X = inter_nonbless, FUN = length)
head(iris.75)Simple feature collection with 6 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 648979.6 ymin: 6861839 xmax: 654353.1 ymax: 6866159
Projected CRS: RGF93 / Lambert-93
CODE_IRIS INSEE_COM geometry nbacc nbaccnb
1 751197316 75119 MULTIPOLYGON (((653970.6 68... 4 2
2 751176716 75117 MULTIPOLYGON (((649189.3 68... 7 3
3 751103703 75110 MULTIPOLYGON (((652767.6 68... 14 5
4 751187104 75118 MULTIPOLYGON (((652827.6 68... 4 2
5 751114314 75111 MULTIPOLYGON (((654272.9 68... 11 4
6 751103707 75110 MULTIPOLYGON (((652960.8 68... 11 5
dplyr : select, group_by et summarize. Ces fonctions fonctionnent également avec les objets sf.
library(dplyr)
com.75 <- iris.75 %>%
group_by(INSEE_COM) %>%
summarize(nbacc = sum(nbacc),
nbaccnb = sum(nbaccnb))
plot(st_geometry(iris.75), col = "ivory3", border = "ivory1")
plot(st_geometry(com.75), col = NA, border = "ivory4", lwd = 2, add = TRUE)Dans cet exercice, nous allons utiliser mapview pour explorer les accidents de la route ayant eu lieu à Paris en 2019.
En complément, nous allons utiliser les données d’OSM permettant de croiser le lieu des accidents avec les routes empruntées.
mapview. Essayez d’utiliser différents paramètres pour customiser votre carte.
map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … du package mapview.
library(mapview)
library(sf)
accidents.2019.paris <- st_read("data/accidents2019_paris.shp")
mapview(accidents.2019.paris, map.types = "OpenStreetMap",
col.regions = "#940000",
label = accidents.2019.paris$Num_Acc,
color = "white", legend = TRUE, layer.name = "Accidents à Paris en 2019",
homebutton = FALSE, lwd = 0.2)Reading layer `accidents2019_paris' from data source `C:\Users\Kim Antunez\Documents\Projets_R\quantilille\exercises\data\accidents2019_paris.shp' using driver `ESRI Shapefile'
Simple feature collection with 11897 features and 7 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 643222.2 ymin: 6857482 xmax: 661451.8 ymax: 6867053
Projected CRS: RGF93 / Lambert-93
Utilisez les polygones de ‘iris.75’ pour extraire la “bounding box” de Paris en projection WGS84.
Récupérez ensuite grâce à osmdata, à l’intérieur de cette bounding box", le fond de carte du périphérique parisien (key = "highway", value = "trunk")
Faites l’intersection entre les accidents et le périphérique, en prenant soin d’ajouter une zone tampon de 50 mètres autour de ce dernier et appeler ce nouvel ensemble de points accidents.2019.paris.periph.
Utiliser sf::st_bbox() et sf::st_transform() pour extraire la bounding box. Le code epsg de WGS84 est 4326.
Utiliser :
osmdata:opq() pour définir la bounding box de la requête osmosmdata:add_osm_feature() pour définir la paire key:value recherchéeosmdata:osmdata_sf() pour récupérer les données osm.
library(osmdata)
#1. bounding box
bb <- iris.75 %>% st_transform(4326) %>% st_bbox()
q <- opq(bbox = bb,timeout = 180)
#2. périphérique
qt <- add_osm_feature (q, key = 'highway',value = 'trunk', value_exact = FALSE)
roads <- c(osmdata_sf(qt))$osm_lines %>% st_transform(st_crs(iris.75))
#3. zone tampon et intersection
accidents.2019.paris.periph <- st_intersection(accidents.2019.paris,
st_intersection(st_geometry(roads),iris.75) %>%
st_buffer(dist = 50) %>%
st_union()
) mapview. Essayez à nouveau d’utiliser différents paramètres pour customiser votre carte.
map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … du package mapview.
library(mapview)
library(sf)
mapview(accidents.2019.paris.periph, map.types = "OpenStreetMap",
col.regions = "#940000",
color = "white", legend = TRUE, layer.name = "Accidents sur le périphérique à Paris en 2019",
homebutton = FALSE, lwd = 0.2)Utiliser la fonction pt_in_grid ci-dessous pour compter le nombre de personnes accidentées dans des cellules de 500m de côté.
mapview pour afficher la grille choroplèthe.
pt_in_grid <- function(feat, adm, cellsize = 1000){
grid <- st_make_grid(x = adm, cellsize = cellsize, what = "polygons")
. <- st_intersects(grid, adm)
grid <- grid[sapply(X = ., FUN = length)>0]
. <- st_intersects(grid, feat)
grid <- st_sf(n = sapply(X = ., FUN = length), grid)
return(grid)
}library(RColorBrewer)
gr <- pt_in_grid(accidents.2019.paris,iris.75,500)
bks = quantile(gr$n)
cols <- brewer.pal(length(bks), "Reds")
mapview(st_as_sf(gr) %>% st_transform(4326),
map.types = "Stamen.TonerLite",
color = "white",
col.regions = cols,
alpha = 0.9,
at = bks,
legend = TRUE,
layer.name = "Nombre d'accidents par</br>carreau de 500 mètres",
homebutton = FALSE, lwd = 0.2)Nous aimerions créer avec le package ggplot2 une carte des arrondissements de Paris qui combine le nombre de personnes accidentées et la part de celles qui n’ont pas été blessées.
Préparation des données :
quantile et RColorBrewer::brewer.pal. Pour la création de la variable typo, vous pouvez utiliser la fonction cut avec les paramètres digit.lab = 2 et include.lowest = TRUE.
library(sf)
# Importer les données
com.75 <- st_read("data/com_75.shp", quiet = TRUE)
# Créer la variable
com.75$part_non_blesses <- 100 * com.75$nbaccnb / com.75$nbacc
# Définir les bks par quantile
bks <- quantile(com.75$part_non_blesses, na.rm = TRUE)
# Définir une palette de couleurs
library(RColorBrewer)
cols <- brewer.pal(length(bks)-1,"Greens")
# For ggplot2 maps - Create a "typo" variable
library(dplyr)
com.75 <- com.75 %>%
mutate(typo = cut(part_non_blesses, breaks = bks,labels = paste0(round(bks[1:(length(bks)-1)])," à ",round(bks[2:length(bks)])),
include.lowest = TRUE))ggplot2, créez une carte qui contient en choroplèthe la variable part_non_blesses et en cercles proportionnels la variable nbacc.
library(ggplot2)
map_ggplot <- ggplot() +
geom_sf(data = com.75, aes(fill = typo), colour = "grey80") +
scale_fill_manual(name = "Part des non-blessés parmi les\naccidentés de la route (en %)",
values = cols) +
geom_sf(data = com.75 %>% st_centroid(),
aes(size = nbacc), fill = "#f5f5f5", color = "grey20", shape = 21,
stroke = 1, alpha = 0.8, show.legend = "point") +
scale_size_area(max_size = 12, name = "Nombre de personnes\n accidentées") +
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(com.75)[c(1,3)],
ylim = st_bbox(com.75)[c(2,4)]) +
theme_minimal() +
theme(panel.background = element_rect(fill = "cornsilk", color = NA),
legend.position = "bottom", plot.background = element_rect(fill = "cornsilk",color=NA)) +
labs(title = "Accidents de la route à Paris en 2019",
caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021") +
guides(size = guide_legend(label.position = "bottom", title.position = "top",
override.aes = list(alpha = 1, color = "#ffffff")),
fill = guide_legend(label.position = "bottom", title.position = "top"))
plot(map_ggplot)reproducibility
sessionInfo()R version 4.0.5 (2021-03-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)
Matrix products: default
locale:
[1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C
[5] LC_TIME=French_France.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggplot2_3.3.3 RColorBrewer_1.1-2 osmdata_0.1.5 mapview_2.9.9
[5] dplyr_1.0.5 sf_0.9-8 knitr_1.33
loaded via a namespace (and not attached):
[1] Rcpp_1.0.6 lubridate_1.7.10 svglite_2.0.0
[4] lattice_0.20-41 leaflet.providers_1.9.0 png_0.1-7
[7] class_7.3-18 assertthat_0.2.1 digest_0.6.27
[10] utf8_1.2.1 R6_2.5.0 leafpop_0.1.0
[13] stats4_4.0.5 evaluate_0.14 e1071_1.7-6
[16] httr_1.4.2 highr_0.9 unilur_0.4.0.9000
[19] pillar_1.6.0 rlang_0.4.11 curl_4.3.1
[22] uuid_0.1-4 raster_3.4-10 jquerylib_0.1.4
[25] rmarkdown_2.7 labeling_0.4.2 webshot_0.5.2
[28] stringr_1.4.0 htmlwidgets_1.5.3 munsell_0.5.0
[31] proxy_0.4-25 compiler_4.0.5 xfun_0.22
[34] pkgconfig_2.0.3 systemfonts_1.0.1 base64enc_0.1-3
[37] htmltools_0.5.1.1 tidyselect_1.1.1 tibble_3.1.1
[40] codetools_0.2-18 fansi_0.4.2 withr_2.4.2
[43] crayon_1.4.1 grid_4.0.5 gtable_0.3.0
[46] jsonlite_1.7.2 satellite_1.0.2 lifecycle_1.0.0
[49] DBI_1.1.1 magrittr_2.0.1 units_0.7-1
[52] scales_1.1.1 KernSmooth_2.23-18 stringi_1.5.3
[55] farver_2.1.0 leaflet_2.0.4.1 sp_1.4-5
[58] xml2_1.3.2 bslib_0.2.4 ellipsis_0.3.2
[61] brew_1.0-6 generics_0.1.0 vctrs_0.3.8
[64] tools_4.0.5 leafem_0.1.6 glue_1.4.2
[67] purrr_0.3.4 crosstalk_1.1.1 yaml_2.2.1
[70] colorspace_2.0-0 rvest_1.0.0 classInt_0.4-3
[73] sass_0.3.1
Iris est un zonage statistique de l’Insee dont l’acronyme signifie « Ilots Regroupés pour l’Information Statistique ». Leur taille est de 2000 habitants par unité.↩︎